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Abstract 

For a system of globally pulse-coupled phase-oscillators, 
we derive conditions for stability of the completely syn- 
chronous state and all possible two-cluster states and ex- 
plain how the different states are naturally connected via 
bifurcations. The coupling is modeled using the phase- 
response-curve (PRC), which measures the sensitivity of 
each oscillator's phase to perturbations. For large systems 
with a PRC, which turns to zero at the spiking threshold, 
we are able to find the parameter regions where multiple 
stable two-cluster states coexist and illustrate this by an 
example. In addition, we explain how a locally unstable 
one-cluster state may form an attractor together will its 
homoclinic connections. This leads to the phenomenon 
of intermittent, asymptotic synchronization with abating 
beats away from the perfect synchrony. 



1 Introduction 

Ensembles of interacting dynamical systems appear as 
mathematical models in many branches of science Gl 0]. 
For example, the study of coupled lasers [j| 0, is 
triggered by technological applications such as high-power 
generation or secure communication 0, HJ . Interacting bi- 
ological units play an important role for the functioning 
of living organisms 0|. Mechanical or electrical coupled 
oscillators have been extensively studied as paradigmatic 
models [Tol 11 1 to study various synchronization phenom- 
ena. In neuroscience, the synchronization of neuron popu- 
lations plays an important role [12L Il3l . Il4j and might even 
lead to pathological effects fl5j |. As a result, there was re- 
cently an increasing effort to control the desynchronization 
of populations of coupled oscillators. In particular, the co- 
ordinated reset stimulation technique [lj, EH proposes to 
establish a cluster-state in the network, in which the oscil- 
lator's phases split into several subgroups. This example 
illustrates the importance of the analysis of cluster forma- 
tion in coupled systems. 



In this paper we investigate the connection between 
the properties of single oscillators and their influence 
on the appearance of clusters in a globally coupled sys- 
tem of such oscillators. We consider the pulse coupling 
(rTI . 18 . [Til . I20I 21 1 approximation, which is widely used 
for modeling of neuron populations. Such an approxima- 
tion is appropriate in the case, when the interaction time 
between the oscillators is much smaller than the character- 
istic period of oscillations. The interaction is mediated by 
the pulses, which are emitted by each of the oscillator after 
reaching some threshold. Although the coupling topology 
is fixed and assumed to be global, i.e. all-to-all and ho- 
mogeneous, the dynamical properties of the individual os- 
cillators will be considered variable (all at the same time) . 
This enables us to describe a bifurcation scenario, which 
links dynamic regimes of stable synchrony with regimes, 
that promote cluster-formation. More specifically, we con- 
sider various sensitivity functions for the individual oscilla- 



tors, called the phase- response-curve [17|,|22j. PRCs have 
been introduced, studied, and computed for many neu- 
ronal models (23, 24, 2^, 22, 26|. They can serve as an ap- 
propriate control parameter determining the properties of 
the oscillators in the network (27| . In our paper, we restrict 
our analysis to PRCs, which turn to zero at the threshold, 
at which the oscillator emits a pulse. This choice is mo- 
tivated by several well known neuron models [l7, 22, 28| 
and means that the uncoupled system at the threshold is 
insensitive to small external perturbations. 

Our work is accomplished by the bifurcation analysis of 
the cluster states and the analysis of the spectral properties 
of the completely synchronous state and eventually emerg- 
ing two-cluster states. Our results reveals how the synchro- 
nization properties of the network depend on the PRC of 
the individual oscillators. We point out how the loss of sta- 
bility of the synchronous solution may give rise to stable 
homoclinic connections of the synchronous solution and in- 
variant two clusters. We derive conditions for the stability 
of these states and, in an exemplary family of unimodal, 
positive PRCs, we illustrate the mechanism by shifting the 
position of their maxima. For some range of the control pa- 
rameter we observe a stable homoclinic connection of the 
synchronous solution, which leads to an apparent synchro- 
nization of the population with eventual beats out of per- 
fect sync at a large time-scale, which are becoming more 
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rare with time. Thus, practically, a synchronized state 
is observed generically in such systems even in the case, 
when the completely synchronous state is locally unstable. 
For another range of parameters, we observe the stabiliza- 
tion of various two-cluster states, which bifurcate from the 
completely synchronous one-cluster state as predicted by 
the analytics. In the paper, the notions "complete syn- 
chronization" and "one-cluster" solution will be used inter- 
changeably as they have the same meaning for our model. 

The structure of the paper is as follows. In Sec. [5J the 
system is introduced. In Sec. [3J we reduce the dynam- 
ics to one-dimensional maps and in Sec. 21 we review the 
stability and bifurcations of the one-cluster state in this 
framework. We also explain how the existence of a stable 
homoclinic set may lead to the phenomenon of intermittent 
asymptotic synchronization. Appearance and stability of 
two-cluster states are studied in Sec. [5] and [51 Numerics 
for some illustrative example is shown in Sec. [7] and some 
technical details are included in the Appendix. 

2 The system 

Networks of weakly coupled oscillators can often be re- 
duced to phase-models, which keep a single phase variable 
if 6 [0,27r] for each oscillator [29|, |22|, |25 1 . In this reduc- 
tion, one naturally encounters a function which measures 
the local sensitivity of an oscillator's phase to small pertur- 
bations - the phase response curve (PRC). If the perturba- 
tions act along only one direction of the state space, e. g. 
the voltage component in many neuron models, the PRC 
can be represented as a scalar function. Further reduction 
is possible when the coupling takes place on a considerably 
smaller time scale than the period of oscillations. In this 
case it is reasonable to approximate the interaction by an 
impact, i.e. by assuming that the interaction is immediate. 
By combining these two reductions, one obtains the model 
of pulse-coupled phase-oscillators 3(1 12, 31, 27 1 , which is 
the subject of our paper. 

We consider the system of N identical phase-oscillators, 
whose dynamics are given as 



tp j (t) = l + -Z(<p j (t)) S(t-t kl ), j = l,...,N. 

(1) 



Here tki are time moments when fc-th oscillator reaches the 
threshold <pk.(tki) = 2-7T, I = 1,2,.... We call these time 
moments also "spikes". At this time, the fc-th oscillator 
sends a spike and all other oscillators j with j ^= fc obtain 
an impact 



C2) 



where Z(<p) is the PRC. At the same time moment, the 
fc-th oscillator resets to (fk = 0. It is assumed that the 
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Figure 1: Examples of different PRCs, (a) Hodgkin- 
Huxley model, (b) Connor model, (adapted from [24|). 



phase response curve Z (if) is smooth inside (0, 27r) , but 
smoothness of derivatives in the endpoints is explicitly not 
assumed, i.e. possibly Z' (0) ^ Z' (2tt) . Moreover, we as- 
sume Z (0) = Z (2tt) = 0, which is frequently met in neu- 
ron models and means that a neuron is insensitive to exter- 
nal stimulation when it is generating or has just generated 
an action potential 17, 22, 28| (see Fig. [TJ. h > is the 
coupling strength. 

An obvious but important property of system (jlj is that 
the order of the phases is invariant if the quotient x/N of 
coupling strength and network size is sufficiently small. 
This means that the oscillators do not overrun each other 
with time. If 



N 
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[0,2tt] Z' {ip)[ 



(3) 



the map ^ becomes strictly monotonous and thus, the 
phase-ordering is preserved during spikes as well. 

Although the dynamical system ([1]) together with the 
resetting condition is defined completely, it is worth to no- 
tice that it is equivalent to a discrete A-dimensional return 
map (tpi, . . . , tpisr) —> K(tpi, . . . , ipn)- This return map can 
be obtained by fixing the position of one oscillator, e.g. 
ifN — 0. Taking into account that the phase ordering is 
preserved, one iteration of the return map corresponds to 
one spike of each oscillator, i.e. N spikes altogether. More 
exact definition of the return map is given in ADpendix l9.il 

Given ([3]), the synchronous solution 



is known to be linearly stable if and only if for all N\ < N 
the following condition holds 



H \ N-Nr 



< 1. 



(4) 



This was proven in [17| by analyzing the linearization of the 
return-map. The terms in are the eigenvalues of this 
linearization at its fixed point ipi = ... = ipjf. Remarkably, 
for the case Z' (0) ^ Z' (2tt) , there might be a loss of 
stability induced by an increasing of the population size. 
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3 Dynamics in invariant cluster 
subspaces 

In this section, we study the dynamics of the return map 
in invariant subspaces of the form 



n 



(fX = ••• = 



ip Nl = 5, 
= <Pn = 



5 G (0,2tt; 



(5) 



In these subspaces, the population is split into two clusters. 
The distance between those clusters is denoted by 6 (or 
2ir — 5), which is the natural coordinate in n^. We can 
reduce the action of return map on IIa^ as follows. Choose 
a particular initial state in iljVn i.e. 



<Pi 



VJVi = 8 an d tp Nl 



+i 



ip N = 0, 



with a particular initial distance 8. Under the dynamics of 
([1]), the next event in time will be the collective burst of 
the oscillators of the first cluster at time t = 2tt — <5. These 
are immediately resetted to (p = 0, i. e. 



ip Nl = 0. 



The impact of the burst on an oscillator of the second 
cluster is obtained by repetitive application of the function 
[i from §2§ - one time for each oscillator in the first cluster. 
This leads to 

<PNi+i = ■■■ = ¥>n = H Nl (27r - 6) , 

where f/ 1 denotes the n-fold su- 
perposition of the function fi, i.e. 
/i"(<5) = /.t(/i(/i(- • • /i(<5) • • • ))). At next, all oscillators 

n 

advance with equal velocity, until the second cluster 
reaches the threshold, that is, 

VNi+l = ■■■ = <£n — 27T. 
Accordingly the first cluster is located at 
2tt - fi Nl (2tt - S) . 



<PN! 



The following burst of the second cluster completes the 
reduction of the return map and after the application of 
the return map, the new distance is given by 

Y Nl (5) = fx N - N ^ (2tt - fj, Nl (2tt -6)), 5e [0, 2tt] . 

(6) 

The one-dimensional map ([6]) determines completely the 
dynamics within the subspace n^, i.e. the dynamics of 
perfect two-clusters ([S]) . Fig. [5] shows typical behaviors 
of these functions and corresponding cobweb-diagrams. A 
fixed point (5* £ (0, 27r) of this map correspond to a two- 
cluster fixed point of the global return map K or, equiva- 
lently, to a periodic two-cluster state of the original system 




s 27t S 2n <L 2n 



Figure 2: Cobweb-diagrams for the one-dimensional maps 
Y/Vj within the two-cluster subspaces 11^. A fixed point 
<5* G (0, 2n) corresponds to a stationary two-cluster state 
of the return map. Chart (a) shows a stable synchronous 
state and an unstable two-cluster state. In (b) one can 
observe a homoclinic connection of the synchronous state, 
resembled by a heteroclinic connection of 6 = and 8 = 2tt 
under , and (c) shows an unstable synchronous state 
together with a two-cluster state, that is stable to per- 
turbations in the subspace T^Nf Possible situations with 
several fixed points for Yn 1 are not shown. In Sec. [SJ we 
explain, how the cases (a) and (b), resp. (b) and (c), are 
connected via the bifurcations of the synchronous solution. 

((T|). In general, the shape of Yjv^J) depends on the clus- 
ter sizes (N\,N — JVi), the PRC Z((p), and the coupling 
strength >c. 

In particular, these subspaces are of interest, since bifur- 
cations along them govern the stability properties of the 
synchronous state, as we show in the next section. 

4 Stability of the synchronous 
state 

The trivial fixed points 6 = and <5 = 2tt of cor- 
respond to the synchronous state. Derivatives of Yn 1 at 
these points can be computed as 

Y' Nl (0) = (l + %Z> (0)) (l + ^'(2-) X 



F^ i (2 7 r)= (1 + _Z'(0)) (l + -Z'(27r)) 



Y^„ Nl (0) 



(7) 



Since the linearized dynamics of the cluster distance 8 for 
small 5 > is governed by 



8^Y^ (0)8, 



(8) 



the derivatives Y^ (0) are multipliers of the map govern- 
ing its local stability at 8 = 0. Hence the linear stability 
conditions of the synchronous state 5 = under Y n for 
n = l,...,N — 1 are |Y^(0)| < 1. These are identical to 
([4]) because the invariant directions ipi = ... = ip^ = 1, 
fNx+i = •-- = ^Pn = are eigenvectors of the return map 
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at the synchronous solution ipi = ... = tpjy = 0. It is ev- 
ident, that the maximal multiplier of the linearized map 
© is either 

y 1 '(o) = (i + ^z'(o)) w_1 (1 + ^(270) 

or 

Yir-i (0) = (l + (0)) 1 (l + %Z' (2tt)) N ~' 

depending on the values of Z'(0) and Z'(2ir). This means, 
that the maximal growth of the cluster distance is observed 
in the case when only one oscillator is separated from the 
rest of the population, i.e. N± = 1 or Ni = N — 1. 

When the number of oscillators N is large, one can ob- 
tain the approximation 

Y Nl (0) ~ exp [x ((1 - p) Z' (0) + pZ' (2tt))] , (9) 

where p = Ni/N denotes the ratio of the oscillator number 
in the first cluster to the total number. This means that 
for large populations, the stability of the synchronous state 
to perturbations in form of two-clusters with N± — pN and 
N — Ni = (1 — p) N elements is governed by the condition 

(1 -p)Z'(0) +pZ'(2tt) < (10) 

and is independent on coupling strength xr. In particular, 
the synchronous solution is linearly stable for large popu- 
lations only if both Z' (2tt) < and Z' (0) < 0. 

5 Two-cluster bifurcations of the 
synchronous state 

In this section, we study generic codimension-1 bifurca- 
tions of the synchronous state, which lead to the emer- 
gence of various two-cluster states. We assume that the 
PRC is varied smoothly in some parameter /3, that is, 

Z {<p) = Zp (ip) . 

In turn, the functions [i (p) = [i (<p, /3) and Yn 1 (5) = 
Y/v-l (5,P) become /3-dependent, too. For the sake of read- 
ability, we will mostly omit this dependency in formulas. 
We still demand that Zp (0) = Zp (2tt) = and that © is 
fulfilled with Z = Zp for all /3. 

As seen in the previous section, the synchronous solution 
is linearly stable in the two-cluster subspace IL^ if and 
only if 

max{Aj Vl (/3),A iV -iv 1 (/3)}<l, (11) 
where 

A„ (p) = Y^(0)=Y N _ n (2ir) 

= (1 + ^(0)) (l+„Z'p(2n)) . 



The loss of stability of the synchronous state in the sub- 
space 11^ takes place if one of the multipliers Ajvi {P) or 
Ajv-Wi (P) becomes bigger than 1 as f3 passes some critical 
value Pi with 

max{AAr 1 (^ x ),Ajv-jVi (fii)} = 1- 

Such bifurcation causes the change of the form of the map- 
ping from those shown in Fig.^a) to Fig.[^b). Note 
that, generally, this loss is only one sided, giving the pos- 
sibility of an emergence of a robust homoclinic connection 
to the synchronous solution as seen in Fig.(2Jb). In a sec- 
ond bifurcation, the stability might be lost completely as 
P reaches a critical value P2 with 

minjAjVi (P2) , Ajv-iVi (P2)} = 1. 

If Pi 7^ /?2, there is either a creation or annihilation of a 
fixed point S t = Y^ 1 (<5„) £ (0, 2ir) in each critical value of 
p. Note that, since Y r [ (2ir) = Y^_ n (0) , the bifurcation 
points P = p c will coincide for N\ = n and Ni = N — n. 
The existence of a fixed point Si of Yn 1 always implies 
that 62 = H Nl {2ir — Si) is a fixed point of Yjv-jVi- There- 
fore, the bifurcations appear as pitchfork bifurcations of 
the return map. 

6 Stability of two-cluster states 

In this section, we will provide explicit conditions for the 
stability of the two-cluster state 

tpi = ... = (p Nl = ip Nl +i = ■•• = Pn = 0. (12) 

These conditions will be given by the formulas (|2T))) . (|2"Tj) . 
and (j2"2")h The practical usefulness of these conditions is 
illustrated later in Sec. [7] by determining the region in 
the parameter space with the coexisting stable two-cluster 
states. Recall, that the two-cluster solution (fl~2|) corre- 
sponds to a fixed point (5* = Yn 1 (5*) . 

Let us investigate the effect of small perturbations to 
p^|) . Because of the original system's symmetry with re- 
spect to index permutations, it suffices to consider per- 
turbations which do not change the phase-order (see Ap- 
pendix [9T] and [l3)- Here, these admissible perturbations 
are of the form 

rj = (7/1, . . . ,r) Nl ,T) Nl+ i, ...,rj N ), 

m> ■■■ > Vn ± > 0, VNt+i > ■ ■ ■ > m = 0. 

To avoid tedious calculations, we do not perform a full 
linearization of the return map and calculation of all mul- 
tipliers. Instead, we will provide estimates, that are sharp 
for N — > 00, for the largest multipliers in the following 
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invariant subspaces of admissible perturbations: 

Vi = {(»yi,.-.,^i-i,0,... J 0) | r?i > • • • > 7?^ = 0} , 

V 2 = {(0, . . . ,0,T) Nl+ i, . . . ,T) N -l,0) 

VN1+1 > ■ ■ ■ > Vn = 0} , 
V3 = {(vi, • • • > Wi , 0, . • • , 0) I 771 = • • • = r) Nl = 6 > 0} 

Here, Vi and V 2 contain inira-cluster perturbations, which 
only affect phases inside one cluster, and V3 contains the 
perturbation of the cluster distance, i. e. the mier-cluster 
perturbations. Since the direct product of these invariant 
subspaces contains all admissible perturbations, the linear 
stability of the full system combines from stability in the 
subspaces. 

Let us first turn to the intra-cluster perturbations in 
Vi and V 2 . We introduce the cluster- widths for the state 
(cp-i, ...,ip N ) as 

Wi=<pi- <pn x and W 2 = PN t +i - Pn- 

In the perfect two-cluster state (fT2")l we have W\ = W2 = 0. 
If this state is perturbed along Vi , the cluster- width of the 
first cluster becomes positive, i. e. W\ = r\\ and W 2 
li. Analogously, perturbations along V 2 result in W\ = 
and W 2 = Vm+i- 

First, let us consider the perturbation within Vi when 
the first cluster with Ni elements is perturbed. The inte- 
stability of one cluster is determined by the change of its 
width during one period, i.e. one application of the re- 
turn map. These changes effectively take place only at two 
events, i. e. at the crossing of the threshold either by the 
first cluster itself or by the second cluster. Denote the max- 
imal possible change of the width of the perturbed cluster 
during its own crossing as A^ 1 (e) , where e > is the 
width before the burst. This is, an initial width e results 
in a new width e + A^ 1 (e) after the burst. Since the other 
cluster is not involved in this process, we can treat it as a 
burst of a single cluster in a system of Ni oscillators. Sec- 
tions H] and [5] showed that, the maximal A^ 1 (e) is realized 
either by the perturbation e = r\\ > 772 = • • • = f/jvi =0 
or by e — r/i = ■ ■ ■ = ?7jvi-i > Wi = 0j depending on the 
values of Z' (0) and Z' (2tt) . 

Similarly, we denote changes of the cluster's width at the 
crossing of the other cluster as A^ - ^ 1 (e) , where e > 
is as before and N — Ni is the number of elements in the 
other, unperturbed cluster. Note that A^~ Nl (e) is the 
same for all intra-cluster perturbations of magnitude e, 
since the phase order is preserved and solely the distance 
between the first and the last phases (fi—tpNi i n the cluster 
matters. 

The maximal width for perturbations of the first cluster 
after one return is then 



Wi (e) = e + Af 1 (e) + A 



N — Ni 



with e = r/i . 

Similarly, the width of the second cluster (within the 
invariant subspace V2) changes as 



W 2 (e) = e + A2 1 (e) + A^" 1 (s + ^ (e)J , 

for e = rjN 1+ i. Stability conditions with respect to all 
possible intra-cluster perturbations (within the subspace 
Vi © V 2 ) are then given by W{ (0) < 1 and W 2 (0) < 1, i. 
e. 



W{(0) = 1+ Af 1 (0)+ A 



L JV-JVi 



(0) 



+ (Af-^)'(0)-(Af)'(0)<l, (13) 

^(0) = l+(A^)'(0)+(Af- Ari )'(0) 

+ (Af-^)'(0).(Af)'(0)<l. (14) 

The stability analysis of the synchronous solution in sec- 
tions [4] and [5l applied to a population of n (n — Ni or 
n = N — Ni) oscillators, implies that A™ (e) satisfies 

1 + (A")' (0) = (l + I min (Z> (0) , Z> (2tt))) ' 



x [1 + — max(Z'(0),Z' (2n)) 

see ([7]). For large populations, this reads 

1 + (A?)' (0) w cxp (^pmax {Z 1 (0) , Z' (2tt))) , (15) 

where p = n/N. This is obtained as in ((9]). The change of 
the width of a perturbed cluster at position tp = 8, that 
is induced by the burst of another unperturbed cluster is 
given as (see Appendix 19.21 for details) 



A 2 ( £ ) = E n ( Z ^ {S + £) ) _ Z ^ {6) ^ ■ (16) 

fe=0 

Hence, 

n 

(A%y(0) = J2^Z'(» k (S))(t* k )'(S)- (17) 



fe=0 



One can approximate (see Appendix 19. 2|) this sum in the 
limit N — > 00 by a simpler expression using the solution 
(r,tp) to the initial value problem 



|f(r^) = xZ(t?(r j¥ >)) ; 
■&(0,<p) = (p. 



(18) 



As argued in the Appendix 19.21 the function (r, ip) is a 
smooth approximation of /j, t ' n ((p) . Using this, we get 



(A?)' (0) 



Z(0(p,6))-Z(6) 
Z(5) 



(19) 
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where p = n/N and 8 is the position of the perturbed clus- 
ter, when the unperturbed cluster crosses the threshold. 

Using (fT5|) . (fl3|) and (fTO|) , we obtain for the two-cluster 
state (fT2| of large populations 



W[ (0) «exp (xpmax (Z' (0) , Z' (2tt))) 

Z(8*) - Z (2k - <d (p,2k - 8*)) 



x 1 



Z(27r-?9(p, 2-K-8*)) 



where p = Ni/N. That is, the condition for the intra- 
stability of the first cluster, (H3J), leads to: 



cxp (xp max (Z' (0) , Z' (2k))) 

( Z (8*) - Z (2tt - d (p, 2tt-S,)) 

X \ + Z(2k-$(j>,2k-8*)) 

and for the second cluster, (|14p reads 

cxp (x (1 - p) max (Z' (0) , Z' (2tt))) 

Z (ti (p,2k - 8*)) - Z (2k - 8*) 



< 1, (20) 



x 1 + 



Z(27T-5*) 



<1. (21) 



The stability of perturbations in V3, i. e. of the cluster 
distance, corresponds to the stability of the cluster position 
<5* as a fixed point of Yat i: i- e. Yj^ (^*) < 1- For JV — >■ 00, 
this condition can be rewritten as 



Z (8*) • Z (■& (p, 2ir-8*)) 
Z (2tt - ?9 (p, 2tt - 8*)) -Z(2k- <5*) 



< 1. 



(22) 



For more details see Appendix [ 

Note that, for x — > 0, i. e. very weak coupling, we can 
simplify ([2H j) -(|2"2 1) to 



(1 -p) Z' (8 )+pmax(Z' (0),Z' (2tt)) < 0, 
pZ' (2k - <5 ) + (1 - p) max (£' (0) , Z' (2k)) < 0, 



(23) 



(24) 



p-Z'(27r-Jo) + (l-p) (5 ) <0, 



(25) 



where <$o = lim^-^o 8* - see Appendix 19.31 These latter 
conditions match exactly the conditions, that were ob- 
tained by Hansel et. al. in [32j for the linear stability 
in the averaged model 



<p i (t) = l + -Y,Z(p i (t)-<p j (t)) 



As in the averaged model, in the pulse-coupled case it 
is possible to encounter multistability in between stable 
cluster-states and the stable synchronous solution. 



6.1 A special case 

Let us now consider the situation 
Z' (0) = Z' (2tt) = 0. 



(26) 



This might seem degenerate, but in neuron models as well 
as in experimental neurophysiology, this property is not 
unusual - review the examples in Fig. Q] and note that 
the derivatives are indeed zero at tp = and tp = 2k. 
We picked a simple example of this kind to illustrate our 
results numerically in the Sec. [7J Given ([26]) . formula ((4]) 
does not supply any information about the stability of the 
synchronous solution, because one finds that for arbitrary 
Ni, 
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and all multipliers are degenerate and equal to 1. Since 
the linear stability analysis does not provide any results, 
one must consider terms of higher order to complete a 
stability analysis for this case. We will go one step in 
this direction by investigating the nonlinear stability of the 
synchronous solution along the invariant split states as §5§ 
and point out that similar mechanisms may participate 
in a destabilization as in the previous section. Consider 
quadratic terms along the split directions. One finds: 

(0) = (N — iVx) ^Z" (0) - N^Z" (2tt) , (27) 

Y& (2k) = (N — N t ) ^Z" (2k) N^Z" (0) 



■Y£- Nl (0) . 



(28) 



Up to the leading terms, the dynamics governing small 
distances between two-clusters is given now by 



1 



(0)«5 2 , 



instead of ©. Thus, we have an analogous situations as 
in Sec. [SJ Fig. [5J The synchronous state is stable under 
conditions 

r JVi (o)<o, y^{2k)>q. 

Bifurcations from stable synchronous state to two- 
cluster-states occur in the same way as before. The condi- 
tions for the stability of two-cluster states (|20|) — (f2T|) from 
Sec. [5] reduce to 



Z(8* 



Z(2n-$(p,2n-8*)) 



< 1 and Z <tfr 27r :?•»<!, 



Z(2k-S*) 



(29) 



and imply (|22p in the present case. Once more, for weak 
coupling, these conditions takes a simpler form, that is 



Z' (So) < and Z' (2ir - 8 ) < 0, 
where 8q = lim^^o 8* . 



(30) 



G 




Figure 3: Family of unimodal PRCs Zp((p) that were used 
in simulations, see 1)331) . 

7 Numerics 

In this section we provide a simple example, where 
branches of initially unstable two-clusters emerge and 
eventually stabilize for increasing /3, as predicted by our 
analytics (see Fig.[S]). Moreover, we will illustrate the sta- 
ble homoclinic structure leading to intermittent synchro- 
nization (see Fig. O . 

In order to detect the appearance of one- or two-cluster 
states, we compute the order parameters 



Ri(t) 



l N 



*¥>j(*) 



and 



N 



(31) 



-Y. 



(32) 



A perfect one-cluster state is characterized by R\ = i?2 = 
1 and a perfect antiphase two-cluster is characterized by 
Ri = and R2 = 1. We present results of simulations for 
k = 0.5, but qualitatively we observe similar behavior for 
a broad range of k > 0. 

Consider the following family of PRCs (see Fig. [3]): 



Z P {if) = 1 - cos fa M) , for /3 e [0, 1] 
(1 



(33) 



X£ (<p) 



2?r 



2/^ 2 



As shown in Fig. 01 we observe two qualitatively dif- 
ferent types of behavior depending on parameter /3. For 
j3 < 0.5, i.e. when the maximum of the PRC is shifted to 
the right (see Fig. [3]), the one-cluster state acts attract- 
ing on most initial states; for /3 > 0.5 the maximum of 
the PRC is shifted to the left and two-cluster states be- 
come attracting. We have chosen initial conditions in a 
vicinity of a two-cluster state in Fig. HJa) and (b) . there- 
fore the initial values of the order parameters are Ri « 
and i?2 ~ 1- Figure 0Jb) shows how the instability of the 
two-cluster state implies desynchronization transient, after 



which the system is attracted to a synchronous one-cluster 
state. Similar behavior occurs for other initial conditions. 
Figure 2Jc) and (d) illustrate the order parameters behav- 
ior for initial conditions close to the splay state (a state, 
where the phases are distributed). The initial values for 
the order parameters in the splay state are close to zero, 
but after a transient, they approach again the same asymp- 
totic values as in (a) and (b). 

A more complicated behavior occurs for the intermedi- 
ate value of the parameter j3 = 0.5, i.e. when the PRC is 
symmetric. In this case, the order parameters Ri(t) and 
i?2(t) do not approach some asymptotic constant values 
but remain periodic in time. As a result, the maximum 
asymptotic values of both R\ and R2 do not coincide with 
the corresponding minimum values. This type of behavior 
is observed for a very small parameter interval of magni- 
tude O {jA around /3 « 0.5. It is worth mentioning, that 
Z0.5 (y) = 1 — cos (ip) was proposed as PRC for oscillators 
near a saddle-node bifurcation on a periodic orbit (24, 22 1. 
Since (3 = 0.5 is a critical value for our model, we advise 
caution for treating this case as exemplary in investiga- 
tions. Figure |4] (e) and (f) summarize the behavior of the 
order parameters for different /3. 

7.1 Cluster-stability 

Since Z' p (0) = Z' p (2tt) = 0, for all j3 e [0, 1] , we have 
Yp,N! (0) = Y^ Nl (2tt) = 1, for all iVj = 1, N — 1. 



This means, expressions (|27[) and (f28[) become relevant to 
stability of the synchronous solutions. In our example, 
they read 

Yp, Nl (0) = 4>,((l- P )/? 2 -p(l-/3) 2 ), 

Y£ Nl (2tt) = 4^ ((1 -p) (1 - /3) 2 - p/3 2 ) , 

where p = Ni/N. For large populations, the synchronous 
solution is unstable for any /3. Nevertheless, numerics show 
an attraction to the synchronous state, which can be ex- 
plained by the existence of a stable homoclinic connection 
of the one-cluster ((33|, see Fig. EJ. 

At j8 ~ 0.5, the symmetric two-cluster state gains sta- 
bility and subsequently, with increasing /?, asymmetric 
cluster-states emerge when Yi' Ni (0) changes sign and gain 
stability when (|29|) is fulfilled. The analytic predictions 
from Sec. |6] match the numerical results. More specifi- 
cally, in Fig. [Ha) one observes pitchfork bifurcations of 
the synchronous state (5 = 0) leading to the appearance 
of two-cluster states. These two-cluster states are initially 
unstable but, with increasing /3, they gain stability creat- 
ing a large set of coexisting stable two-clusters. Figure[6][b) 
shows the stable clusters versus f3. Solid line corresponds 
to the analytically predicted stability domain (|29|) and the 
dots are observed clusters from numerical calculations. 
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Figure 4: Asymptotic behavior of the order parameters. 
Charts (a) and (b) show Ri(t) and i?2(i) for a trajectory 
starting in a vicinity of the two-cluster state. The mid- 
dle panel (c) and (d) belong to a trajectories starting in 
a vicinity of the splay state. Left charts (a) and (c) cor- 
respond to the parameter value 8 = 0.7, where the two- 
cluster state is attracting and (b) and (d) to 8 = 0.3, where 
one-cluster state is attracting. In (e) and (f): dependence 
of the asymptotic values for the order parameter R\ and 
i?2 on 8. For the most values of 8, except 8 = 0.5, the 
order parameters tend to some constant value, when ini- 
tialized near the splay state or the symmetric two-cluster 
state. 




5.000 



5.000 



Figure 5: Homoclinic connection of the one-cluster. In 
(a): Synchronous solution (i.e. one-cluster state) as 
a saddle point in the phase space with a homoclinic 
loop. In (b) and (c): Width of the cluster A(t) = 
maxi<i .j<N {\fi{t) — Vj^W as a function of time. Chart 
(a) shows the behavior along the orbit started at an initial 
condition close to the splay state (far from the one-cluster) 
and (b) shows the behavior along the orbit started close to 
a split-state with homoclinic connection as in Fig. [2] (b). 
The growing intervals between peaks in (a) indicate the 
existence of a stable homoclinic orbit. 




Figure 6: Stability and existence of two-cluster states, (a) 
A cascade of pitchfork bifurcations gives birth to fixed 
points of Yn 1 and Yjv-iVi- Solid lines denote stable two- 
cluster stationary states and dashed - unstable. The lines 
are shown only for selected values of p = Ni/N. Markers in 
Fig. (b) shows which two-clusters are stable in dependence 
on 8 and p (observed numerically). The value p = 0.5 
corresponds to the symmetric cluster and p ^ 0.5 to non- 
symmetric clusters. The shaded region in (b) indicates 
the analytic prediction for the stability of the two-cluster 
states by 
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8 Discussion 

In this paper, we have considered the following question: 
How the individual properties of oscillators in a globally 
coupled system influence the formations of one- and two- 
cluster states? Working in the framework of pulse-coupled 
oscillators, the adjustable parameter for our study was cho- 
sen to be the PRC. 

Starting from the results of Goel and Ermentrout [ItJ 
about the stability of one-cluster state, we have extended 
these results to the case when the first derivatives of the 
PRC at the threshold equal to zero. We have derived a one- 
dimensional map that describes the dynamics in invariant 
subspaces corresponding to various two-cluster states. By 
means of this map we identified pitchfork bifurcations of 
the one-cluster state, which lead to the emergence of peri- 
odic two-cluster states. For large populations, we provide 
explicit conditions for the linear stability of these states. 
The analysis of the one-dimensional map also reveals that a 
homoclinic connection to the synchronous one-cluster state 
is generic in the considered class of systems. Moreover, nu- 
merical analysis shows that a set of homoclinic connections 
can be globally attracting. 

For the numerical illustration, we have chosen systems 
with positive, unimodal PRCs, which turn to zero at the 
threshold together with its first derivatives. This corre- 
sponds to the excitatory coupling. By varying the shape 
of this PRCs, we observe how the initially globally at- 
tracting one-cluster undergoes a sequence of bifurcations 
in which two-cluster states emerge and later on stabilize. 
This leads to the coexistence of multiple stable two-cluster 
states. The latter phenomenon resembles qualitatively the 
Eckhaus phenomenon 34. 35|. 
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9 Appendix 

9.1 Formal definition of the return map 

In this section, we derive the discrete return map for the 
system (|T|). Let us introduce the state-spaces 

x j = {[<Pi: -j Vn) € [0, 2?r] I (p [1+j] > ... > <P[ N+ j] = 0} , 

where [n] = n mod N. In particular 

X N = X = {(<pi,...,tp N ) € [0,2tt] I (p! > ... ><p N = 0}. 



Since all oscillators are taken to be identical and homo- 
geneously coupled, the system is symmetric with respect 
to permutations of indexes (©Ar-symmetry). Therefore, 
without loss of generality, let us consider initial states in 
Xq. The return map K : Xq h-> Xq is composed of N maps 
kj : Xj-i i y Xj, j = 1, N, defined as 

kj {(p 1 ,...,(p j ,...,(p N ) := 

<p 1 +2ir- (pj + —Z {<pi +2w-ipj),... 



, ip N + 2n - tpj + —Z (tp N + 2-k 



<Pi. 



where the j — th The map kj describes the firing of the 
j-th oscillator, assuming it has maximal phase. Under 
the assumption that the phase order is preserved, i.e. @ 
holds, the return map is given by 



K = A; 7v o ■•• o ki : Xq — > X± 



X 



N 



x . 



Note that each kj is smooth. Linearizations of K capture 
perturbations, which keep the presumed phase-ordering. 
Because of the system's symmetry, it suffices to consider 
these for a determination of linear stability (see also (l7j). 

9.2 Approximation of repetitive firing 

Assume that ^ holds. Then, the phase order is preserved. 
Consider a perturbed cluster, ipi > • • • > tpx 1 , with N\ = 
p ■ N and 

<pi = 6 + e, tpNi = 5. 

The change of the width of this cluster during the threshold 
crossing of another unperturbed cluster 

VJVi+i = ■ ■ ■ = ifN = 2tt 

with N2 = (1 — p) N oscillators is given by the change of 
the distance between the leading oscillator ipi and the last 
oscillator y?jy, . The new position ipf of a single oscillator 
ipj of the perturbed cluster is given as 

<pf = m W2 (vj) = KKK- ■ ■ Kfj) •••))) 



N 2 

N 2 -l 



fe=0 



Thus, the change of width of the perturbed cluster is 



ft - ftr, 



n 2 -i \ 

k=0 ) 



N 2 -l 



-( 5+ E ^(m*(*)) 



JV2-1 



fc=0 
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In this way, we arrive at (|16p . 

Let us introduce a function, # : [0, 1] x [0, 2ir] — > [0, 27r] , 
that approximates the iterations of fj, as follows: 



■d (r, tp) ~ p fe (tp) , for r 



N' 



Then, 



r + -,tp)-# (r, tp) « M r - 7V+1 (V) - (*>) 



AT 



Z(0(r,¥>)), 



or equivalently, 



fl(r + ^,y) -fl(r,y>) 
1/AT 



xZ(i?(r, V )). 



Therefore, for large N , ■& (r, tp) is a solution to the initial 

d± 

dtp 



value problem (JTHJ. Further, ^p(r, solves the linear 



system 
5 dti 



ch) 



drdtp ^,tp)=xZ'('d(r,t P )) — (r,tp), 



(0,^ = 1. 



Therefore, 

0^ (r, <p) = cxp I J xZ ($(s,tp))ds 
Taking into account the property (|18|l . we obtain 



f^ )=cxp (Z 



*(r,*0 z / (g) \ Z (fl( r , y )) 



(34) 



Now, we can simplify the expression (|17[) by means of 
For this, we substitute ^ k (tp) by #(-^,<p) and the sum 
over fc by an integral over r in (|17|) . As a result, we obtain 



AT 2 -1 
k=0 

xZ' (0(r, 5)) — (r, 5) dr 
dip 

Z(d(l-p,S))-Z(5) 

Z(5) 

Thus, for large populations, (|17p leads to (|19|). 

The inter-cluster-stability is given by the stability of 5* 
as a fixed point of Yp.N, that is, |l^.jy (<5*)| < 1- In terms 
of 1? (r, v?) , we have Y p . N (S) i? (1 - p, 2tt - 1? (p, 2tt - <5)) 
and fixed points <5* = Y p . n (<5») satisfy J* w 
1? (1 -p, 2tt - 1? (p, 2tt - d*)) . Thereby one obtains ([2"2"j) . 



9.3 Weak coupling in large populations 

Assume that for sufficiently small x, there exist fixed 
points 5 (x) = Y p . N (S (x) , x) , with p ■ N e {1, . . . , N} . 
Further, assume that the limit So = lim^^o 8 ( K ) exists. 
Actually, it is a generic property for any point S g (0, 2n) , 
that fulfills (1 — p) Z (So) = pZ (2tt — So) , that there exists 
a family of fixed points of Y p .n with So = lim >r _ s .o & (x) ■ 

For small x, the linearization of (|20| - ((22|) in x = 
gives the stability conditions (|23| -([25 |) . For example, the 
linearization of (l20l) is 



0>^-( exp(xpmax(Z'(0),Z' (2tt))) 

<7X \ 

1+ Z{2*-<Hp,2*-6.)) " 



In the following, we add x explicitly as an argument 
where needed. Note that the two-cluster fixed points 5* = 
5 (x) are x-dependent as well as the function $ (r, tp) = 
$ (r, tp, x) from the previous section and Y p .n (tp) = 
Yp.N (tp, x) . Using (fl8| . we obtain gjj (0, tp, x) as a solu- 
tion to 

9 dd 



dr d> 



(r,tp,x) = Z($(r,tp,x)) 



+ xZ' (d (r, tp, x)) — (r, tp, x) , 
ox 



(0,<p,x) = 0. 



This can be solved explicitly as 

d-d 

— (r, tp,x)=r ■ Z (i? (r, tp, x)) . 



(36) 



Using t?(r,<p,0) = tp, (J34J) and ([36]). the condition (J35J) in 
x = becomes 



> pmax(Z'(0),Z'(27r)) 

d (Z(S(x))-Z(2-K-§(p,2-K-S(x),x)) 
dx 



pmax(Z' (0),Z'(2ir)) 



Z(2ir-$(p, 2tt-S(x),x)) 

pZ (2tt - So) Z' (So) 



Z(So) 



(37) 



Differentiation of the fixed point equation S (x) = 
Y p . N (S (x) , x) w (1 - p, 2tt - t? (p, 2tt - 5 (x) , x) , x) in 
x = yields 



w (l-p)Z(<Jo)-pZ(27r-5o). 



(38) 



Inserting this into ([3T)) gives ([2^|) . Similarly, one obtains 
from dH]) and J25J from 
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